Métodos Estadísticos

• La modelización mediante métodos estadísticos permite

  • Encontrar patrones

  • Interpretar datos

tibble(x = runif(n=1e7, min = 0, max = 20)) %>% pull(x) %>% mean()

set.seed(9915)
tibble(x = runif(n=25, min = 0, max = 20)) %>% pull(x) %>% mean()

Procedimiento Bootstrap (cont.)

distribDF %>% ggplot(aes(x=medias)) + geom_density()

# Error estándar de la media
mediaSE <- distribDF %>% pull(medias) %>% sd()

# Intervalos de confianza al 95% para la media
mediaIC <- distribDF %>% pull(medias) %>% quantile(c(0.025, 0.975))

Aprendizaje supervisado (cont.)

• Dos formas de especificar la relación \(f\)

  • Modelo paramétrico: supone un forma de \(f\) que depende de parámetros desconocidos, p.e., lineal \(f(x) =\beta_0 + \beta_1 x_1 + \dots + \beta_k x_k\)

  • Modelo no paramétrico: ajustar \(f\) a los datos sin supuestos funcionales

    • más flexibilidad y mejor ajuste, pero más difícil de estimar e interpretar

• La capacidad predictiva de un modelo mejora si incluimos más variables explicativas (modelo más flexible)

• Dependiendo de la naturaleza de la variable dependiente, tenemos:

  1.- Problemas de Regresión: variable de respuesta cuantitativa (toma valores numéricos)

  2.- Problemas de Clasificación: variable de respuesta cualitativa (toma valores en una de \(C\) categorías o clases)

Regresión Lineal: Estimación

• Objetivo: estimar los coeficientes poblacionales desconocidos a partir de una muestra

• El error de estimación o residuo es \(\small \hat{e}_i = y_i - \hat{y}_i\), donde la predicción a partir del modelo estimado es \(\small \hat{y}_i=\hat{\beta}_0+\hat{\beta}_1 X_1+ \dots + \hat{\beta}_k X_p\)

• Los coeficientes estimados son los que minimizan la Suma Cuadrática de Residuos: la suma total de distancias entre los datos observados y predichos

\[ \small SCR=\sum_{i=1}^{n} \hat{e}_i^2= \sum_{i=1}^{n} ( y_i - \hat{y}_i)^2 \]

• Esto equivale a las condiciones derivadas de suponer \(\small E(\varepsilon|X)=0\)

Regresión Lineal: Estimación (cont.)

• Los parámetros del modelo \(f(.)\) estimado, en este caso los \(\small \beta_j\) son estadísticos con variabilidad muestral, medida por su error estándar \(\small se(\widehat{\beta}_j)\)

• Se usan para construir intervalos de confianza y para contrastar hipótesis sobre los parámetros poblacionales, p.e., significatividad (\(\small \beta_1=0\))

• En R, usamos la función lm()

descuento <- read.csv("data/descuento.csv") %>% 
              mutate(zona = parse_factor(zona))
modelo <- lm(data = descuento, ventas ~ renta + descuento)
summary(modelo)

• La predicción \(\widehat{y}\) también está sujeta a incertidumbre por la estimación: se puede calcular su error estándar e intervalos de confianzas

ventas_pred <- predict(modelo, type = "response")
cbind(descuento$ventas, ventas_pred) %>% head()

Regresión Lineal: superando la linealidad

• Se pueden incluir transformaciones no lineales de las variables del modelo

• La interpretación del cambio de un regresor sobre \(\small Y\) es diferente

lm(data=descuento, log(ventas) ~ renta + I(renta^2) + descuento) %>% summary()

• Otra forma de permitir no linealidades es discretizando variables continuas: permite efectos “escalón” diferentes para distintos valores

lm(data=descuento , ventas ~ cut_width(renta, 500) + descuento) %>% summary()

• Se incluyen indicadores binarios para cada clase excepto uno
  • la constante recoge el valor medio de \(\small Y\) para ese grupo de referencia
  • cada coeficiente recoge el efecto adicional de su grupo sobre \(\small Y\)

Regresión Lineal: superando la linealidad (cont.)

• También podemos incluir interacciones entre variables: el efecto de un regresor dependerá de otro regresor

lm(data=descuento, ventas ~ renta*descuento)  %>% summary()

• Principio jerárquico: al incluir una interacción siempre deben incluirse los factores principales (NO sólo renta:descuento)

• Cuando interactuamos un regresor continuo y uno binario, permitimos que la pendiente del primero sea diferente para cada grupo

lm(data=descuento , ventas ~ (renta + zona)*descuento) %>% summary()

–>

“Problemas” del Modelo de Regresión Lineal

• No linealidad: incluir transformaciones no lineales

• Correlación de los errores: afecta a los errores estándar, no la estimación

  • usar errores estándar robustos o modelizar la dinámica

• Heterocedasticidad: ídem, usar errores estándar robustos

  • los gráficos de los residuos frente a un regresor o valores predichos: ¿heterocedasticidad o no linealidad?

• Outliers en la variable de respuesta o en los regresores

• Colinearidad: indica que no es posible separar el efecto de cada regresor: eliminar alguno o recombinarlos

• No normalidad: TCL, Bootstrap,…

• El único supuesto realmente importante es \(\small E[\varepsilon|X]=0\)

Regresión Logística

• La regresión lineal puede usarse respuestas binarias (no más de dos categorías), aunque genera predicciones fuera del rango \(\small [0,1]\)

• Solución: aplicar al índice lineal una transformación \(\small F(z)\in[0,1]\)

• De manera que \(\small \Pr(Y=1|X)= p(x)= \Lambda( \beta_0 + \beta_1 x_1 + \ldots + \beta_k x_k)\)

Regresión Logística (cont.)

• Los coeficientes NO se interpretan como cambios en la probabilidad ante cambios unitarios en un regresor (efecto marginal sobre la probabilidad)

• PERO su signo (y significatividad) son los mismos que los del efecto marginal

• Como NO tiene sentido minimizar la SCR, el objetivo es maximizar la probabilidad (verosimilitud) de observar los unos y ceros en los datos

• La regresión logística pertenece a la familia de modelos lineales generalizados (GLM, en inglés)

• Se pueden incluir como variables explicativas tanto variables cuantitativas como cualitativas, e incurrir en sesgo por omisión de variables

glm(data = Default, default ~ student, family = "binomial" ) %>% summary()
glm(data = Default, default ~ student + balance, family = "binomial" ) %>% summary()

Regresión Logística: Predicciones

• El objeto de R de glm() contiene valores predichos, que son probabilidades de \(\small Y=1\)

logit <- glm(data = Default, default ~ balance*student, family = "binomial" ) 
cbind(Default$default, logit$fitted)

• También se puede predecir usando una muestra distinta de la usada para estimar o con valores concretos de los regresores

logit <- glm(data = Default, default ~ balance, family = "binomial" )
predict(logit, newdata = tibble(balance=c(0,100)), type="response")

Regresión logística con más de dos clases

• La regresión logística se puede generalizar a situaciones con múltiples clases (modelos multinomiales) con un índice lineal para cada clase \[ \small \Pr(Y=c|X)=\frac{e^{\beta_{0c}+\beta_{1c}X_1+\dots+\beta_{kc}X_k}}{\sum_{l=1}^{C}e^{\beta_{0l}+\beta_{1l}X_1+\dots+\beta_{kl}X_k}} \]

• La librería glmnet() permite la estimación de estos modelos

library(glmnet)
x <- model.matrix(Species ~ Sepal.Length + Sepal.Width, data = iris)
mod.glmnet <- glmnet(x = x, y = iris$Species, family = "multinomial", 
                     lambda = 0, type.multinomial = "grouped")
coef(mod.glmnet) 
predict(mod.glmnet, newx=x, type = "response")  # probabilidad de cada clase
predict(mod.glmnet, newx=x, type = "class")     # clase

Error de predicción

• Un modelo es mejor si sus predicciones se ajusten mejor a las observaciones

• El error de predicción es \(y - \widehat{y} = f(X) - \widehat{f}(X) + \varepsilon\)

  • \(f - \widehat{f}\) = error reducible (eligiendo modelo)

  • \(\varepsilon\) = error irreducible (variables no observadas)

• La función de pérdida (o coste) evalúa cómo valoramos las desviaciones

Métricas de error de predicción (cuantitativa)

• Mean Square Error (Error Cuadrático Medio): \(\small MSE(y,\widehat{y})={\frac{1}{n}\sum_{i=1}^{n}\left(y-\widehat{y}\right)^2}\)

  • penaliza grandes desviaciones

  • \(\small R^2\) y \(\small R^2\)-ajustado son variantes del MSE, y solo sirven para comparar modelos con la misma variable dependiente.

• Root Mean Square Error: \(\small RMSE(y,\widehat{y})=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(y-\widehat{y}\right)^2}\)

  • mismas unidades que \(\small y\)

• Mean Absolute Error: \(\small MAE(y,\widehat{y})=\frac{1}{n}\sum_{i=1}^{n}\left|y-\widehat{y}\right|\)

• Otras medidas basadas en distintas funciones de pérdida, la verosimilitud del modelo (\(\small AIC\), \(\small BIC\)), etc

Seleccionar el mejor modelo

• ¿Podemos predecir el número de visitantes en función de la temperatura?

library(mosaicData)
data(RailTrail)
RailTrail %>% ggplot(aes(x = hightemp, y = volume)) + 
  geom_point() +  geom_smooth()

• ¿Cuál es el mejor modelo?

  • \(\small volume = \beta_0 + \beta_1 hightemp + \varepsilon\)

  • \(\small volume = \beta_0 + \beta_1 hightemp + \beta_2 hightemp^2 + \varepsilon\)

  • …

  • \(\small volume = \beta_0 + \beta_1 hightemp + \dots + \beta_{22} hightemp^{22} + \varepsilon\)

RailTrail %>% ggplot(aes(x = hightemp, y = volume)) + 
  geom_point() + 
  geom_smooth(method = 'lm', formula = y ~ poly(x,22) ) + 
  coord_cartesian(ylim = c(100,750))

Muestras de entrenamiento y de prueba

• Las métricas de error (ej., \(\small MSE\)) suelen calcularse para predicciones de los mismos datos usados para ajustar/estimar el modelo (in-sample prediction)

  • Esta muestra se denomina muestra de entrenamiento (training sample)

• PERO queremos saber qué tal se predicen casos nuevos (out-sample prediction)

• Usar las métricas en muestras de entrenamiento implica problemas de “overfitting”: sobreajuste a las características de la muestra concreta

  • Un modelo menos flexible podría tener menor error de predicción con casos nuevos

• Debemos calcular las métricas de error con observaciones que el modelo NO ha usado antes: muestra de prueba (test sample)

“Overfitting”

• Siempre que aumenta la flexibilidad, el MSE

  • disminuye en la muestra de entrenamiento

  • tiene forma de U en la muestra de prueba

• Nota: el MSE en entrenamiento es siempre menor que en prueba

MSE en la muestra de prueba

\[ \small E\left[\left(y-\widehat{f}(x)\right)^2\right] = E\left[\left(f(x)-\widehat{f}(x) + \varepsilon + E\left[\widehat{f}(x)\right]-E\left[\widehat{f}(x)\right] \right)^2\right] = \]

\[ \small =\underbrace{\left[E\left(\widehat{f}(x)\right)-f(x)\right]^2}_{(1)} + \underbrace{E\left(\left[\widehat{f}(x)-E\left(\widehat{f}(x)\right)\right]^2\right)}_{(2)}+Var(\varepsilon) \]

• \(\small (1)=\left[Sesgo\left(\widehat{f}(x)\right)\right]^2\): error por supuestos erróneos en \(f\)

  • ajuste insuficiente (underfit) al perder relaciones relevantes entre \(X\) e \(Y\)

• \(\small (2)=Var\left(\widehat{f}(x)\right)\): sensibilidad a fluctuaciones en el entrenamiento

  • si el algoritmo modela puro ruido en entrenamiento, ajustará bien allí, pero predecirá mal casos nuevos (overfit)

“Trade-off” Varianza–Sesgo

• El sesgo se reduce y la varianza aumenta con la complejidad del modelo \(\Rightarrow\) encontrar un método (ej., flexibilidad) para el que ambos sean bajos

• NO es posible minimizar simultáneamente ambas fuentes de error: memorización (en entrenamiento) vs. generalización de resultados

Medir el Error en la Clasificación

• Los modelos de clasificación NO predicen directamente la categoría, sino la probabilidad de que una observación pertenezca a cada categoría

• La clase predicha será aquella con mayor probabilidad. En el caso binario, implicar superar el umbral de 0.5 (se deben probar varios valores)

censo <- read_csv("data/census.csv") %>%
  mutate(income = parse_factor(income))
logit <- glm(income ~ capital_gain, data = censo, 
             family = "binomial")
prob.predict <- predict(logit, type = "response")

umbral <- 0.5
cat.predict  <- if_else(prob.predict > umbral, 1, 0) 
cbind(censo$income, cat.predict, prob.predict) %>% head(10)

• Como no tiene sentido diferencia de clases (variables categóricas), NO se pueden calcular medidas como el MSE y otros relacionados

Matriz de Confusión

• Tabular categorías observadas frente a las categorías predichas

table(cat.predict, censo$income)

• Existen varias medidas derivadas de la matriz de confusión

Métricas con la matriz de confusión

• Tasa de observaciones correctamente clasificadas (exactitud o accuracy)

\[ \scriptsize ACCUR=\frac{VP+VN}{VP+FP+VN+FN} \]

• No es informativo cuando algunas clases son infrecuentes (datos desequilibrados)

  • si hay poco fraude/enfermos (ej., 5%), predecir que nunca hay fraude implica \(\scriptsize ACCUR=95\%\), PERO NO detecta fraude/enfermedad

• El estadístico Kappa (\(\small \kappa\)) es una medida similar, pero que ajusta por lo se esperaría solo por azar (corrigiendo en parte el desequilibrio entre clases).

Métricas con la matriz de confusión (cont.)

• La tasa de verdaderos positivos o sensibilidad (recall) es el porcentaje de verdaderos positivos sobre el total de positivos observados \[ \scriptsize TVP=SENSIT=\frac{VP}{VP+FN} \]

  • ej., tasa de fraude/enfermos existentes que se detectan correctamente

• La tasa de verdaderos negativos o especificidad es el porcentaje de verdaderos negativos sobre el total de negativos observados \[ \scriptsize TVN=ESPECIF=\frac{VN}{VN+FP} \]

  • ej., tasa de “otras” opciones que se clasifican correctamente

  • Tasa de falsos positivos: \(\scriptsize TFP = 1 - TVN = 1 - ESPECIF\)

Métricas con la matriz de confusión (y 3)

• La exactitud equilibrada (Balanced Accuracy) es una media de la sensibilidad y de la especificidad

• La precisión o valor de predicción positivo es la cantidad de verdaderos positivos sobre el total de positivos predichos

\[ \scriptsize PREC=\frac{VP}{VP+FP} \]

• La familia de medidas \(\small F_{\beta}\) es una ratio de la importancia ponderada de la sensibilidad y de la precisión: \(\scriptsize F_{\beta}=\frac{(1+\beta)^2 \times SENSIT \times PREC}{\beta^2 \times SENSIT + PREC}\)

  • Para \(\scriptsize \beta<1\), se da menos importancia a la sensibilidad: los falsos positivos se consideran más costosos

  • Para \(\scriptsize \beta>1\), los falsos negativos son más costosos y para \(\scriptsize \beta=1\) son igualmente costosos

Curva ROC (“Receiver Operating Characteristic”)

• Representa TVP (eje y) frente a TFP (eje x) en diferentes umbrales : reducir el umbral clasifica más elementos como positivos (verdaderos y falsos)

• La curva ROC informa del grado de separabilidad: dado un nivel de TFP, el clasificador es mejor cuanto mayor sea TVP

AUC (“area under the curve”)

• La AUC es el área bajo la curva ROC: ofrece una medida agregada de rendimiento entre 0 (todas las clasificaciones incorrectas) y 1 (todas correctas)

• Resume la curva ROC y permite comparar curvas que se cruzan

Extensiones. Métricas adicionales

• Con más de dos clases, se realiza un análisis AUC-ROC para cada categoría (frente a las demás) y se promedian (ej., ponderando por casos en cada clase)

• Con clases desequilibradas, se puede preferir en lugar de la ROC un gráfico de precisión frente sensibilidad (precision-recall) y su correspondiente AUC (PR-AUC)

• Existen múltiples funciones de pérdida (o coste de clasificación) posibles.

  • Las relacionadas con la curva de ganancia consideran el coste de alcanzar un cierto nivel de sensibilidad

  • Otras se basan en la función de verosimilud o la entropía como medidas de pérdida (ej. mean log loss)

Evaluación de Modelos: entrenamiento y prueba

• Para minimizar problemas de underfit y, sobre todo, de overfit, DEBEMOS dividir aleatoriamente el conjunto de datos en dos partes:

• Entrenamiento (80-90%): datos sobre los que se construye/estima el modelo

• Prueba(20-10%): se usa el modelo construido para predecir y se evalúa con datos no vistos antes

• ¿Por qué renunciar a parte de los datos si sabemos que un tamaño muestral grande es importante? Evaluar correctamente un modelo lo es mucho más

• La estimación del error en prueba puede ser volátil dependiendo de las observaciones incluidas en cada grupo

Evaluación de Modelos: Validación cruzada

• Para evitar que los datos sean sensibles a una partición concreta, se usa validación cruzada (cross-validation o rotation estimation)

• Se repite varias veces y de forma ordenada el proceso de remuestreo para la partición en grupos de entrenamiento y prueba (similar a bootstrap)

• Permite utilizar todas las observaciones de la muestra, tanto para estimar como para evaluar el modelo (aunque no a la vez)

Validación cruzada de K bloques

• Se divide, aleatoriamente y ex-ante, la muestra en K subconjuntos (normalmente, 5 o 10)

• Un subconjunto se usa como prueba y el K-1 restantes como entrenamiento

• Se repite el proceso durante k iteraciones, con cada posible subconjunto de datos de prueba.

• Se obtiene una métrica de error en cada iteración; se promedian para obtener un único resultado de evaluación

• Es el tipo más habitual de validación cruzada

Validación cruzada aleatoria (RCV) y LOOCV

• RCV: en cada iteración se realiza la particion aleatoria (con reemplazamiento) entre entrenamiento y prueba

• Las observaciones pueden “repetir” como prueba

• LOOCV (leave one out CV): solo una observación se usa como prueba en cada iteración y el resto como entrenamiento

• Se realizan \(n\) iteraciones; se calcula una media sobre \(n\) resultados